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ABSTRACT 

As many as 5 ice giants — Neptune-mass planets composed of ~90% ice and rock and ~10% 
hydrogen — are thought to form at heliocentric distances of ~10-25 AU on closely packed orbits spaced 
^5 Hill radii apart. Such oligarchies are ultimately unstable. Once the parent disk of planetesimals 
is sufficiently depleted, oligarchs perturb one another onto crossing orbits. We explore both the onset 
and the outcome of the instability through numerical integrations, including dynamical friction cool- 
ing of planets by a planetesimal disk whose properties are held fixed. To trigger instability and the 
ejection of the first ice giant in systems having an original surface density in oligarchs of S ~ 1 g/cm 2 , 
the disk surface density a must fall below ~0.1g/cm 2 . Ejections are predominantly by Jupiter and 
occur within ^10 7 yr. To eject more than 1 oligarch requires a < 0.03 g/cm 2 . For certain choices of a 
and initial semi-major axes of planets, systems starting with up to 4 oligarchs in addition to Jupiter 
and Saturn can readily yield solar-system-like outcomes in which 2 surviving ice giants lie inside 30 AU 
and have their orbits circularized by dynamical friction. Our findings support the idea that planetary 
systems begin in more crowded and compact configurations, like those of shear-dominated oligarchies. 
In contrast to previous studies, we identify a < 0.1S as the regime relevant for understanding the 
evolution of the outer solar system, and we encourage future studies to concentrate on this regime 
while relaxing our assumption of a fixed planetesimal disk. Whether evidence of the instability can 
be found in Kuiper belt objects (KBOs) is unclear, since in none of our simulations do marauding 
oligarchs excite as large a proportion of KBOs having inclinations > 20° as is observed. 

Subject headings: celestial mechanics — Kuiper belt — planets and satellites : formation — solar system 
: formation 



1. INTRODUCTION 

Without gravitational focussing, in situ coagulation of 
Uranus and Neptune takes too long to complete. In a 
minimum-mass disk at heliocentric distances of 20-30 
AU, timescales to assemble the ice giants exceed the age 
of the solar system by 2 orders of magnitude, if growth 
is unfocussed (e.g., Goldreich, Lithwick, & Sari 2004, 
hereafter GLS04). 5 N-body coagulation simulations that 
do not damp relative velocities between planetesimals, 
either by dynamical friction, inelastic collisions, or gas 
drag, fail to form Uranus and Neptune (Levison & Stew- 
art 2001; see also Lissauer et al. 1995). The ice giants 
contain 10-20% hydrogen by mass, a fraction so large 
that such gas must originate from the solar nebula. The 
outer planets must therefore form within a few x 10 7 yr, 
before all of the nebular hydrogen photoevaporates (Shu, 
Johnstone, & Hollenbach 1993; Matsuyama, Johnstone, 
& Hartmann 2003). 

One way to alleviate (but not necessarily eliminate) 
the timescale problem is to form Uranus and Neptune 
closer to the Sun, where material densities and collision 
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rates are greater. Thommes, Levison, & Duncan (1999, 
2002) explore a scenario in which the two planets form 
at distances of 5-10 AU, between the cores of Jupiter 
and Saturn. Once the gas giant cores amass their en- 
velopes, they scatter the ice giants outward onto eccentric 
orbits. These orbits subsequently circularize by dynami- 
cal friction with planetesimals at 15-30 AU. Tsiganis et 
al. (2005) propose an alternative history in which Uranus 
and Neptune accrete at 12 and 17 AU, are thrown out- 
ward by Jupiter and Saturn, and have their orbits circu- 
larized by dynamical friction. According to their story, 
the outward scattering of ice giants is triggered by hav- 
ing Jupiter and Saturn divergently migrate across their 
mutual 2:1 resonance. 

Another approach to solving the timescale problem is 
to consider how gravitational focussing can be amplified. 
GLS04 adopt this route by appealing to a massive disk 
of sub-km-sized planetesimals, similar to those produced 
by coagulation simulations set in the outer solar system 
(Kenyon & Luu 1999). The disk envisioned by GLS04 has 
a mass several times the minimum-mass value in conden- 
sates so that the "isolation mass" — the mass to which a 
edu protoplanet grows by consuming all material within its 
annulus of influence — equals Neptune's mass. The small 
bodies comprising the disk collide so frequently that their 
velocity dispersion damps to values less than the Hill ve- 
locity of a Neptune-mass planet. Accretion rates then 
enjoy maximal enhancement by gravitational focussing, 
and proto-Neptune can accrete the last half of its mass 
in ^10 5 yr (see eqn. 105 of GLS04). Gas drag supplies 
another means to damp planetesimal velocity dispersions 
(Rafikov 2004; Chambers 2006; see also appendix A of 
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GLS04). 

Strongly focussed, in situ assembly of planets from a 
dynamically cold disk carries, however, a potential em- 
barrassment of riches: The disk can spawn more ice 
giants than the solar system's current allotment of 2 
(Uranus and Neptune). We estimate that about 5 iso- 
lation masses or "oligarchs," each having the mass of 
Neptune, can form between 15 and 25 AU (see eqn. Q] 
below). These planets comprise a "shear-dominated oli- 
garchy," so-called because the encounter velocities be- 
tween planets and planetesimals are given by their min- 
imum values set by Keplerian shear. Initially, the oli- 
garchs' nested orbits would be stabilized by dynamical 
friction with the disk. GLS04 suggest that excess oli- 
garchs would be purged from the outer solar system by an 
eventual dynamical instability. According to their order- 
of-magnitude analysis, this "velocity instability" occurs 
once the mass of the disk becomes less than the mass 
in oligarchs, whereupon dynamical friction ceases to sta- 
bilize the system against mutual gravitational stirring 
(a.k.a. "viscous stirring"). In the ensuing chaos, sev- 
eral oligarchs would be ejected, either by other oligarchs 
or by Jupiter or Saturn, possibly leaving two survivors 
whose orbits could circularize by dynamical friction at 
15-30 AU. 

Despite their disparate perspectives on the timescale 
problem and different motivations, the scenarios of 
Thommes et al. (1999), Tsiganis et al. (2005), and GLS04 
share quite a few features. In their simplest forms, each 
theory starts with a more crowded configuration for so- 
lar system planets than is observed today; each is char- 
acterized by an intermediate period of dynamical chaos 
during which Uranus and Neptune execute highly eccen- 
tric orbits; and each invokes final regularization of ice 
giant orbits by dynamical friction with an ambient disk. 
Thommes et al. (2002) and Chiang et al. (2006, here- 
after C06) point out that these violent histories might 
be encoded in Kuiper belt objects (KBOs). In particu- 
lar, so-called scattered KBOs possess large eccentricities, 
inclinations, and perihelion distances which might reflect 
gravitational stirring by marauding ice giants. 

The notion that planets originate in compact and 
crowded configurations is bolstered by the study of extra- 
solar systems as well. To explain the striking preponder- 
ance of large orbital eccentricities observed among extra- 
solar giant planets (Butler et al. 2006), multiple planets 
each having on the order of a Jupiter mass are imag- 
ined to have once resided on orbits sufficiently close that 
the planets scatter one another onto elliptical trajecto- 
ries (Marzari & Weidenschilling 2002; Ford, Rasio, & Yu 
2003; Ford, Lystad, & Rasio 2005). 

GLS04 outline a possible formation history for Uranus 
and Neptune in a packed oligarchy, and C06 expand upon 
its consequences for the Kuiper belt, by making many 
simplifying assumptions and order-of-magnitude approx- 
imations. In this paper, we test some of their ideas by 
numerical simulations. In particular, we seek answers to 
the following questions: 

1. For a shear-dominated oligarchy containing more 
than two Neptune-mass oligarchs beyond Saturn's 
orbit, what is the critical value of the disk surface 
density below which the velocity instability occurs? 

2. What is the likelihood that the instability will re- 



sult in the survival of two oligarchs whose final or- 
bits resemble those of Uranus and Neptune? 

3. To what degree is the Kuiper belt dynamically ex- 
cited by velocity-unstable oligarchs? 

Our methods are described and tested in <J2] That 
section contains empirical determinations of how rapidly 
5 oligarchs viscously stir one another, with and with- 
out the gas giants Jupiter and Saturn. Comparisons are 
made with analytic theory. Results of hundreds of sim- 
ulations designed to provide statistical answers to the 
above questions are presented in Ej3j We summarize and 
offer an outlook in $4j 

2. METHOD AND TESTS 

To guide the reader, we provide a condensed descrip- 
tion of our method in £| 2 . H Details are elaborated upon 
in ^2.2H2.3I 

2.1. Overview 

We simulate the final stages of oligarchy by numerically 
integrating the trajectories of 5 closely packed Neptune- 
mass oligarchs, together with those of 2 gas giants resem- 
bling Jupiter and Saturn. Oligarchs and gas giants are 
referred to as planets. We employ the hybrid orbit in- 
tegrator MERCURY6 (Chambers 1999), which combines 
a conventional Bulirsch-Stoer integrator to handle close 
encounters between planets, with the fast symplectic al- 
gorithm invented by Wisdom & Holman (1991). Viscous 
stirring of an oligarch by other planets is simulated as ac- 
curately as the orbit integrator solves the gravitational 
equations of motion. Case studies of viscous stirring are 
described in £12.21 

To model dynamical friction between a planet and the 
surrounding disk of planetesimals, we introduce a pertur- 
bative force on each planet. The force damps the com- 
ponent of the planet's velocity that differs from the local 
disk (circular) velocity. For simplicity, we take the disk 
to have a constant surface density between an inner and 
an outer radius. Disk parameters are held fixed. In our 
simple scheme, the planets respond to the disk through 
dynamical friction, but the disk does not respond to the 
planets. The details of the perturbation force are pro- 
vided in H2.3I The validity of our fixed disk approach is 
briefly considered in £ )4.2.21 

Finally, to investigate how oligarchs might excite the 
Kuiper belt, an ensemble of test particles is included in 
a subset of the simulations. We use the terms "test par- 
ticle" and "Kuiper belt object (KBO)" interchangeably. 
These test particles are intended to represent large KBOs 
like those observed today, having sizes on the order of 100 
km. This size is small enough that we can neglect dynam- 
ical friction between KBOs and the disk, yet also large 
enough that we can ignore damping of KBOs' velocities 
by physical collisions with the disk (C06). Thus, in our 
simulations, KBOs (test particles) feel directly only the 
gravity of the Sun and of the planets. 

In the sub-sections below, we explore separately the 
processes of viscous stirring f §2.2p and dynamical fric- 
tion ( §2.3| . in isolation from one another. We present 
full-fledged simulations, in which the two processes are 
combined, in Sj3j 
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2.2. Viscous Stirring 

We first study how multiple oligarchs gravitationally 
stir one another. For this sub-section, we ignore dynam- 
ical friction with the disk, but include the gas giants, 
Jupiter and Saturn. We compare our findings to those 
of GLS04. The results of this sub-section will be applied 
in fJ3]to understanding the threshold conditions required 
for velocity instability. 

2.2.1. Initial Conditions 

We consider iV iig = 5 oligarchs, each having the mass 
of Neptune (fj, = m/rriQ = mx/m® = 5.1 x 10~ 5 ). The 
oligarchs are initially spaced 5 Hill radii apart in semi- 
major axis (a); that is, the difference between semi-major 
axes of nearest-neighboring oligarchs is 

a j+1 - aj = 2.5(/x/3) 1/3 a 3 + 2.5(^/3) 1/3 %+i , (1) 

where j ranges from 1 to iV ii g . The coefficient of 2.5 is 
inspired by numerical studies by Greenberg et al. (1991) 
of the width of a protoplanet's feeding zone, for the case 
where Keplerian shear dominates the relative velocity be- 
tween a protoplanet and a planetesimal (see their eqn. 9; 
see also Ida & Makino 1993). While the coefficient of 2.5 
is the standard value for our study, the coefficient in re- 
ality can be somewhat different, depending on the accre- 
tion and migration histories of the planets. We explore 
some consequences of varying the coefficient in £j2.2.3l and 

BH 

We assume the semi-major axis of the innermost oli- 
garch a\ = 15 AU. Then according to eqn. ([1]), semi- 
major axes of the next four oligarchs equal 17.1, 19.4, 
22.1, and 25.1 AU. Initial eccentricities, and initial incli- 
nations relative to an arbitrary reference (x-y) plane, are 
such that ej = smij = 10 -4 . All orbital elements in this 
paper are osculating and referred to a barycentric coor- 
dinate system. That is, when computing the osculating 
Kepler elements for a given body, the position and ve- 
locity of the body are referred to the system barycenter 
(calculated using all massive bodies), while the central 
mass of the assumed Kepler potential equals the mass of 
the Sun plus that of the given body (alone). For each oli- 
garch, the initial longitude of ascending node, longitude 
of pericenter, and mean anomaly are randomly generated 
from uniform distributions between and 2-7T. 

Jupiter and Saturn are assigned initial masses and 
semi-major axes equal to their current values: fij = 
mj/m e — 9.5 x 1CT 4 , = ms/m = 2.9 x 1CT 4 , 
aj = 5.18 AU, and as = 9.54 AU. Initial eccentricities 
equal ej = es = 0.05, and inclinations relative to the ref- 
erence plane are such that sinij = sinig = 0.01. Orbital 
longitudes are randomly generated, just as they are for 
the oligarchs. 

These initial conditions, particularly our choices for 
ai = 15 AU and N Q n g = 5, are somewhat arbitrary. They 
are intended to represent qualitatively the final stages 
of shear-dominated oligarchy in the outer solar system. 
We will adjust starting parameters (e.g., ai, AT iig) in 
later sections to achieve simulation outcomes in better 
agreement with observed properties of the solar system. 

A total of iV rea ] = 200 orbital integrations ("re- 
alizations" ) are performed with the hybrid integrator 
MERCURY6 (Chambers 1999), each characterized by a 



unique set of starting longitudes and each lasting 10 7 yr. 
The timestep for the symplectic integrator is set to 130 
days. Timesteps for the conventional Bulirsch-Stoer in- 
tegrator are as short as necessary to achieve an accuracy 
parameter of 10~ 10 . The changeover distance that sep- 
arates the symplectic regime from the close encounter 
regime is set to Ar cr i t = 3 Hill radii. 

A "collision" between two bodies occurs when their 
mutual separation becomes less than the sum of their 
physical radii. The physical radius of each oligarch is 
computed using Neptune's bulk density of 1.6g/cm 3 . 
Physical radii for Jupiter and Saturn are computed using 
densities of 1.3 and 0.7g/cm 3 , respectively. We assume 
that bodies that collide merge completely 

An "ejection" occurs when an oligarch's distance from 
the Sun exceeds 10000 AU and when its total kinetic 
plus potential energy (evaluated in barycentric coordi- 
nates with the potential energy set to zero at infinity) 
becomes positive. Ejected planets are dropped from the 
simulation. 

2.2.2. Results: Outcomes After t = 10 7 yr 

At t — 10 7 yr, the outcomes for all A a n = A^ ea ] x 
Aoiig = 200 x 5 = 1000 oligarchs divide into the follow- 
ing mutually exclusive categories, in order of decreasing 
frequency of incidence: 

1. Ejection but no collision (463) 

2. No collision and no ejection (439) 

3. Collision with another oligarch but no ejection (42; 
i.e., 21 oligarchs remain, each with mass twice that 
of an original oligarch) 

4. Collision with Sun only (19) 

5. Collision with Jupiter only (14) 

6. Collision with another oligarch, and subsequent 
ejection or subsequent collision with Jupiter or Sat- 
urn (13) 

7. Collision with Saturn only (10) 

The dominant outcome is ejection. In 50% of the re- 
alizations (i.e., 100 out of -/Vreai — 200), at least one 
oligarch is ejected before t = 1.6 x 10 6 yr. By t — 
3.2 x 10 6 yr, 85% of all realizations experience a first ejec- 
tion. All but 5 out of 200 realizations experience at least 
one ejection of an oligarch by t = 10 7 yr. 

Jupiter and Saturn are responsible for the preponder- 
ance of ejections. When we repeat the experiment with 
Jupiter and Saturn omitted, outcomes at t = 10 yr are 
as follows: 870 out of 1000 oligarchs experience neither 
an ejection nor a collision, 118 collide with another oli- 
garch (so that 59 remain), and 12 collide with two other 
oligarchs (so that 4 remain). No ejection is observed to 
occur by t = 10 7 yr when the gas giants are absent. 

These outcomes are consistent with timescale esti- 
mates by GLS04. Neglecting Jupiter and Saturn, GLS04 
predict that the oligarch ejection timescale is ^10 9 yr (see 
their eqn. 114). This is consistent with our finding that 
no ejection occurs within t = 10 7 yr in the absence of 
gas giants. GLS04 mention the possibility that Jupiter 
and Saturn hasten ejections. We confirm this possibil- 
ity. When gas giants are present, we find the ejection 
timescale is ^10 6 yr. 
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the mass surface density of oligarchs. Conjunctions 
with these oligarchs occur over the synodic period 



1 10 100 1000 10" 10 5 10" 10 7 
Time (years) 

Fig. 1. — Viscous stirring of oligarchs, in the absence of dy- 
namical friction, for an oligarch separation of 5i?H- We show the 
median (solid curve), 85 th percentile (long dashed curve), and 15 th 
percentile (dotted curve) of the oligarchs' eccentricity distribution 
versus time, based on N-body integrations with no dynamical fric- 
tion. In panel (a), each integration includes Jupiter, Saturn, and 
five oligarchs. In panel (b), only the five oligarchs are integrated. 
Arrows in panel (a) mark the times when 15%, 50%, and 85% of 
iV rea ] = 200 realizations (integrations) experience the first ejection 
of an oligarc h. Ver tical lines separate various phases of evolution 
discussed in i|2,2.3l 

2.2.3. Results: Eccentricity and Inclination Growth 
("Stirring Curves") 

Fig. Q] tracks the median eccentricity, eso (i) , of all oli- 
garchs. The sample from which the median is drawn al- 
ways contains iV a n = 1000 objects, regardless of whether 
oligarchs collide or are ejected. When computing the 
median, we adopt the following rules: ejected oligarchs 
have their eccentricities set equal to 1 (but remain part 
of the sample); an oligarch that collides with either the 
Sun, Jupiter, or Saturn has its eccentricity set equal to 
1; and oligarchs that collide with other oligarchs are still 
counted as separate objects and are each assigned an ec- 
centricity equal to the current eccentricity of the merged 
body. Experiments with alternative sets of rules pro- 
duced no qualitative changes to our results. 

The resultant "stirring curves" of Fig. [T] exhibit a vari- 
ety of behaviors. We first discuss the case when Jupiter 
and Saturn are omitted from the integrations (Fig. [TJa) . 
As annotated in Fig. Q}), we distinguish four phases of 
viscous stirring: 

1. Distant Conjunctions: At early times t < 500 yr, 
planetary orbits do not cross and eso grows roughly 
linearly with time. A linear dependence is expected 
from viscous stirring by distant conjunctions, i.e., 
conjunctions between oligarchs that are not nearest 
neighbors. To derive the t 1 scaling, we estimate us- 
ing the impulse approximation that a conjunction 
between two oligarchs separated by distance x < a 
imparts eccentricities on the order of 



(2) 



to both bodies, provided they have eccentricities 
less than Ae prior to conjunction. For a given 
oligarch, a total of N ~ Y.ax/m oligarchs all re- 
side about the same distance x away, where X is 
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as roughly observed in Fig. Q] As time elapses, 
ever closer neighbors at smaller x drive the stir- 
ring. This reasoning matches that given by GLS04 
in their treatment of viscous stirring in the shear- 
dominated regime; their eqn. (33) is identical in 
form to our eqn. ((4]). The t 1 scaling is also derived 
by Collins & Sari (2006) and Collins, Schlichting, 
& Sari (2007). These latter studies concentrate on 
the limit iV ii g ^> 1. 

2. Conjunctions with Nearest Neighbors: At inter- 
mediate times 500 < t(yr) < 1.5 x 10 4 , plane- 
tary orbits remain non-crossing but the eccentric- 
ity distribution hardly changes. Since the syn- 
odic period between nearest neighboring oligarchs 
is t syn ~ 500 yr, a given oligarch during this phase 
experiences repeated conjunctions with its nearest 
neighbor. Such repeated close encounters might be 
expected to produce chaotic motion and to cause 
eccentricities to random walk, in which case eso oc 
i 1 / 2 . That this scaling is not observed implies that 
our 5 oligarchs do not behave in a strongly chaotic 
manner despite their close spacing. Indeed, we ob- 
serve in our simulations that the epicyclic phases 
(true anomalies) of a given oligarch at successive 
conjunctions with its nearest neighbor do not vary 
completely randomly. Perturbations from conjunc- 
tions with a nearest neighbor apparently tend to 
cancel out during this second phase. 

3. Onset of Orbit Crossing: The cancellations charac- 
terizing the preceding phase are not perfect, how- 
ever. Eventually, from t ~ 1.5 x 10 4 yr to t ~ 
5 x 10 4 yr, eccentricities surge as oligarchs start 
crossing orbits. The median eccentricity eso sur- 
passes the orbit-crossing value, e w 0.06, during 
this third phase. 

Our finding that orbits cross in a few x 10 4 yr is 
consistent with numerical experiments by Cham- 
bers, Wetherill, & Boss (1996), who measure times 
required for close encounters to occur in initially 
circular, co-planar, multi-planet systems as a func- 
tion of planet mass and orbital spacing. For refer- 
ence, our spacing of 5i?n corresponds to 4 "mutual 
Hill sphere radii" as defined by those authors. 

4. Orbit Crossing: At late times t > 5 x 10 4 yr, oli- 
garchs routinely cross orbits and we observe eso c< 
^0.25 -y^ e can re p r oduce this scaling using the fol- 
lowing particle-in-a-box argument. An oligarch's 
random velocity v at time t is determined largely 
by its closest encounter with another oligarch up 
until that time. We call the impact parameter 
characterizing this closest encounter 6 m i n . From ki- 
netic theory, nb^ n vt ~ 1, where n ~ YiQ,/mv is the 
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number density of oligarchs and we have assumed 
that the random velocity distribution of oligarchs 
is isotropic. It follows that b m ; n oc 1/t 1 / 2 and 
v ~ (Gm/b 2 ^) 1 / 2 oc t 1 / 4 . This scaling agrees with 
that of eqn. (49) of GLS04 and that of eqn. (14) of 
C06. 

When we restore Jupiter and Saturn to the integrations 
(Fig.QJi), we can still discern the four phases enumerated 
above. However, compared to the case without giants, 
some phase boundaries are shifted to earlier times, and 
eso rises more quickly during some phases. Phase 1 tran- 
sitions to phase 2a at t ss 30 yr; at this time, all oligarchs 
have undergone their first conjunctions with Jupiter and 
Saturn. Phase 2a transitions to phase 2b at t w 500 yr; as 
in the case without giants, this transition marks the time 
when every oligarch has experienced about one conjunc- 
tion with its nearest neighboring oligarch. During phase 
2b, we witness the same remarkable near-constancy of 
e 50 « 0.01. Finally, during phase 4 at t > 3 x 10 4 yr, 
when oligarchs are on crossing orbits, e§o oc t 8 . Such 
growth outpaces that observed in the absence of the gas 
giants. 

What about oligarch inclinations? In simulations with- 
out gas giants, we observe that the median inclination 150 
remains fairly constant at the initial value of 10 until 
phase 3. As orbits cross, sin 150 surges up to ~0.05, and 
thereafter grows as t 025 during phase 4, just as eso does. 
By t — 10 7 yr, i 50 « 10°. When Jupiter and Saturn are 
included, sini 50 oc t - 28 during phase 4. By t — 10 7 yr, 
i 50 w 10°, as was the case without the gas giants. The 
modest growth of oligarchs' inclinations will limit the de- 
gree to which inclinations of KBOs are stirred f £)3.3.2[) . 

Not all of the different phases of viscous stirring that 
we observe are anticipated from the study of GLS04, 
which documents only the t 1 scaling characterizing shear- 
dominated oligarchy (phase 1) and the t 1 / 4 scaling char- 
acterizing the super-Hill, orbit-crossing regime (phase 4, 
no giants). Their analysis misses the intermediate phase 
2 of slow-to-no growth just prior to orbit crossing, and 
the significant roles that Jupiter and Saturn play in ac- 
celerating viscous stirring during phase 4 (t - 38 vs. t 1 / 4 ). 
That differences exist is not too surprising, as their anal- 
ysis is rooted in the large iV ii g S> 1 limit, whereas for our 
system, iV iig = 5. More importantly, we take nearest- 
neighboring oligarchs to be separated by 5 Hill sphere 
radii, as dictated by the extent of an oligarch's feeding 
zone in a shear- dominated disk (Greenberg et al. 1991), 
whereas the order-of-magnitude equations of GLS04 gov- 
erning shear-dominated oligarchy assume the separation 
is closer to ~1 Hill sphere radius. In this regard, we 
present in Figure [2] viscous stirring curves for cases where 
the oligarch separation is 3i?H and 7i?n (corresponding 
to coefficients in eqn. (fTJ of 1.5 and 3.5, respectively). 
Without gas giants, for a separation of 3i?n, phases 2 
and 3 disappear, leaving only phases 1 and 4 as origi- 
nally envisioned by GLS04. The time to orbit crossing 
varies from ~300 yr to ~2 x 10 6 yr as the spacing changes 
from 3 to 7 Hill radii. This extreme sensitivity to spac- 
ing was also found by Chambers et al. (1996). Including 
gas giants, however, reduces this sensitivity, as Fig. [2^ 
shows. 

The actual oligarchic spacing might only be determined 
by careful numerical simulations of accretion and orbital 
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Fig. 2. — Viscous stirring of oligarchs, in the absence of dy- 
namical friction, for assumed oligarchic spacings of 3-Rh an d 7-Rh- 
Data are generated and presented in the same way as for Fig. [T] 
but only median eccentricities are displayed. Without the gas gi- 
ants, the same 4 phases characterizing viscous stirring for 5Hh 
(Fig. [TJ are evident for 7-Rh- A spacing of 3/?h produces systems 
sufficiently chaotic that phase 2b is indiscernible. Adding the gas 
giants, however, reduces the sensitivity of the results to the oli- 
garchic spacing. 

migration. We adopt in this paper a standard value of 
5i?n, identical to that assumed by GLS04, and motivated 
by studies of shear-dominated accretion by Greenberg et 
al. (1991). Shorter spacings seem less attractive insofar 
as they will produce smaller isolation masses for a given 
disk surface density. 

The intermediate phase of slow-to-no growth of eccen- 
tricity that characterizes oligarchic spacings > 5i?n will 
prove important in determining the threshold disk sur- 
face density below which dynamical friction cooling can- 
not balance viscous heating, i.e., the threshold surface 
density for the velocity instability f £|3.2p . 

2.3. Dynamical Friction 

Oligarchs grow from a disk of planctcsimals. Those 
planetesimals that are not accreted exert dynamical fric- 
tion on oligarchs. The conditions for velocity instability, 
and the ease with which survivors of the instability return 
to low-eccentricity, low-inclination orbits, depend on the 
strength of dynamical friction. We describe how we im- 
plement dynamical friction in our simulations in §2.3.11 
present a test case in ^2.3.21 and show that our imple- 
mentation is compatible with the formulae of GLS04 in 

2.3.1. Prescription 

Consider a planet having an eccentricity and an in- 
clination much greater than those of disk planetesimals. 
Dynamical friction reduces the planet's random (pecu- 
liar) velocity: the difference v = vv between the orbital 
velocity of the planet and that of the mean disk flow. 
From Binney & Tremaine (1987, their eqn. 7-17), 



dv 
~dt 



2irG 2 mp 



m(l + A 2 )£, 



(5) 



where p is the local mass density of the disk, m is the 
mass of the planet, G is the gravitational constant, and 
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A = 



2 (sin 2 «disk)^, 2 



)^circ) 



Gi 



(6) 



is the Coulomb parameter appropriate for dynamical fric- 
tion in a Keplerian disk (Stewart & Ida 2000). For 6 max , 
the maximum impact parameter between the planet and 
a disk planctesimal, we adapt the expression of Stewart 
& Ida (2000; see the discussion following their eqn. 2-17): 



Vax = Ru + r ((sin 2 i disk ) 



• 2 a V 2 
sin i) 



(7) 



where r is the instantaneous distance between the planet 
and the system barycenter, i?n = W^) 1 ^ 3 ?" is the Hill 
sphere radius, and (sin 2 idisk) 1 ^ 2 <C 1 is the inclination 
dispersion of disk planetesimals, held constant for each 
simulation (more on its precise value later). The term 
2(sin 2 «disk) u ci rc in eq n - @ approximates the square of 
the velocity dispersion of disk planetesimals, where w c i rc 
is the local mean disk speed. We take u c i rc to equal 
the speed that the planet would have on a circular orbit 
about the Sun. 

Usually it is assumed in writing eqn. ([5]) that A > 1. 
We do not make this assumption. In fact, we use 
eqns. ©-© regardless of the magnitude of A. When 
A <C 1, dynamical friction is in the shear-dominated 
regime. In £ j2.3.3[ we justify our universal application 
of ©-([71) by showing that these equations correctly re- 
duce to forms appropriate to the shear-dominated case 
when A -C 1 . 

We implement dynamical friction as follows. We are 
interested in the case where the oligarchs are so dynami- 
cally excited that each plunges through a vertically thin, 
dynamically cold disk of planetesimals twice per orbit. 
Specifically, we assume that the time a planet spends 
immersed in the disk, At ~ h/\v z \, where h is the full 
vertical thickness of the disk and \v z | is the vertical com- 
ponent of v at the moment of disk crossing, is short com- 
pared to the orbital period, i or b = 27r/fL Equivalently, 
sini 3> (sin 2 idisk) 1 ^ 2 - At every disk crossing, a planet 
receives a specific impulse of 



Av 



dv 

~dl' 



- At « j — m (! + A 2 )pAt v 



■ ln(l + A 2 )- — -v. 



(8) 



where a is the disk surface density (height-integrated p) . 
At every timestep of the integration, we check whether 
the planet crosses through the disk, which is fixed to lie 
in the x-y plane. At moments of disk crossing, we apply 
a kick according to eqn. ©: we increment the velocity of 
the planet by Av but do not change the planet's position. 
We compute the difference velocity v by subtracting the 
barycentric velocity of the planet from v c - lrc (j), where </> 
is the unit vector in the azimuthal direction. The kick 
is applied in the subroutine MDT_HY.FOR in the MER- 
CURY6 code, after the positions are advanced but before 
the velocities are updated for the second time by the in- 
teraction Hamiltonian. 

For all our simulations, we fix (sin 2 idisk) 1//2 = 10 -3 , 
a value sufficiently small that sini ^s> (sin 2 idisk) 1 ^ 2 for 
all but a tiny fraction of the time. In other words, the 



strength of dynamical friction in our simulations depends 
much more strongly on the planet's random speed v than 
on the much smaller random speeds of planetesimals (see 
eqns. [Br©. Planetesimals can maintain low velocity dis- 
persions by inelastic collisions or by gas drag. 

Our scheme for dynamical friction damps orbital in- 
clinations relative to the x-y (disk) plane. The inclina- 
tion may become so small that At <x 1/ sin i exceeds 
t or b, at which point the planet is immersed within the 
disk and our impulse approximation breaks down. To 
account for this possibility, we arbitrarily set At = 
vam{h/\v z \, 0.025/O). Our softening prescription slows 
but does not stop the damping of inclination and eccen- 
tricity for sin i < 0.004. The softening might represent 
slight misalignments between the planet's orbital plane 
and the disk midplane, which in reality will be warped. 
We have verified that our principal findings, described in 
$3[ do not depend sensitively upon the details of this 
prescription. While the precise values of the inclina- 
tions that we compute are clearly not very meaningful, 
we expect that our results are still qualitatively correct, 
i.e., the code correctly identifies when mutual inclina- 
tions between planetary orbits are large (> lrad) and 
small (<C 1 rad). 

The main virtue of our prescription for dynamical fric- 
tion is its simplicity. We need only specify the disk sur- 
face density cr(r), not the volumetric density p(r, z), and 
we need only apply dynamical friction at disk crossings. 
The main shortcoming of our prescription is that it does 
not account for the response of the disk to the plan- 
ets. The clearing of gaps and the generation of time- 
dependent, non-axisymmetric structures (see, e.g., Gol- 
dreich & Tremaine 1982) will alter the gravitational back- 
reaction of the disk onto embedded planets in ways that 
could be significant but that we (and GLS04) do not cap- 
ture. For ways to model the response of the planetesimal 
disk more realistically, see Lithwick & Chiang (2007) and 
Levison & Morbidelli (2007). 



2.3.2. Test with Single Planet 

We test our prescription for dynamical friction in the 
case of a single planet interacting with a disk. The disk 
has constant surface density a = lg/cm 2 . A Neptune- 
mass oligarch is placed on an orbit having initial semi- 
major axis ajnit = 20 AU and initial eccentricity and 
inclination such that e^t = sinzj n i t = 0.3. Two cases are 
considered, one where the initial argument of perihelion 
Wi n it = and another where Wmit = 7r/2. The evolution 
of e(t) and i(t) depends on Winit mod n. 

Fig. [3] displays the resultant evolution. The planet's 
eccentricity and inclination both drop precipitously to- 
ward zero after a time on the order of 10 6 yr; the exact 
time varies by a factor of 3 between our choices for u init . 
The semi-major axis can increase or decrease. It changes 
by 3-13%, on the order of but less than the starting ec- 
centricity. 

We check our numerical results by comparing them to 
the following approximate analytic solution. Since the 
kick Av is applied twice per orbit, we write 



dv 
~di 



2Av 

^orb 
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Fig. 3. — Test of our prescription for dynamical friction. We 
calculate the orbital evolution of a single Neptune-mass planet in- 
teracting with a disk of surface density <r = 1 g/cm 2 . Solid curves 
denote the evolution when the planet's initial argument of peri- 
center oJinit = (so that the orbit intersects the disk at pericenter 
and apoocenter), while dashed curves correspond to u>i n i t = tt/2 
(so that the orbit intersects the disk at quadrature). Dotted curves 
represent our analytic solution JTTJ which assumes a constant semi- 
major axis a and Coulomb parameter A. 



-47rln(l + A 2 )- 



torb V 3 \V Z \ 



(9) 



We set e = sinz and make the following approximations: 
v = (ep + iz)fla, v — \^2efla, and \v z \ = ifla. Here p and 
z are unit vectors that lie parallel and perpendicular to 
the disk, respectively. Then ^ simplifies to 



de _ ln(l + A 2 ) Ga m 1 
dt y/2 n,a m© e 3 



(10) 



with an analogous equation for i. For fixed A and a, 
eqn. (fT0|) integrates to 



41n(l + A 2 ) to Ga 
y/2 wq fla 



1/4 



(11) 



The eccentricity (equivalently, inclination) vanishes in a 
finite time 



t 



V2 



vanish 



41n(l + A 2 ) to Ga 



(12) 



We overlay eqn. (fTT|) in Fig. [3J taking as representative 
values a = a init = 20 AU and ln(l + A 2 ) = ln(l + A 2 nit ) = 
13. The analytic solution lies between the two numerical 
solutions. We consider the agreement acceptable. 



2.3.3. Connecting Our Prescription to GLS04 

Our equations can be re-cast into the same forms as 
those of GLS04, under the assumption e = i. We start 
with eqn. ([9]) and substitute \v z \ = t OI h = 2n/tt, 

and m = (47r/3)p p i?p, where p p and R p are the internal 
density and physical radius of the planet: 



ldv 
v dt 



-2 3/2 ln(l + A 2 ) afl 



4-kG 2 Pp RI 



3v 4 



(13) 



We next recognize that 47rG 2 j o 2 i?p/3 = 3Vg Scp /T67r, 
where v CS c,p is the escape velocity from the surface of 
the planet. Then eqn. (fT3|) simplifies to 



1 dv 
v dt 



3\/2 



ln(l + A 2 



i aQ 



(^esc, 
V 



(14) 



When ln(l + A 2 ) is a constant of order unity, eqn. (I14| 
matches the form of eqn. (45) of GLS04, evaluated using 
the first line of their eqn. (46), with their planetesimal 
random velocity u replaced by v (since v > u; see their 
section 5.5, end of first paragraph). This formula de- 
scribes dynamical friction in the dispersion-dominated 
regime, where v exceeds the Hill velocity «h = ^-Rh- 

On the other hand, it is possible for A <C 1. This 
happens, according to when v ~ ^/2ifl <g; 

\J Gm/Rft ~ «h (terms proportional to (sin 2 idisk) 1 ^ 2 = 
10~ 3 are negligible). In this shear-dominated regime, 



ln(l + A 2 ) w A 2 , and eqn. ([14j) reduces to 



(15) 



(16) 



ldv _ 3 an /i? H 
v dt tt\/2 Pp-Rp V -Rp 
3 an 1 
ny/2 PpRp a 2 ' 

where we have defined, following GLS04, a = R p /Rn- 
Eqn. p6p matches, to within a numerical constant, 
eqn. (45) of GLS04, evaluated using the second line of 
their eqn. (46). 

We conclude that our treatment of dynamical friction 
is compatible with that of GLS04. 

3. RESULTS 

We present the results of simulations that combine vis- 
cous stirring due to multiple oligarchs, with dynamical 
friction due to a planetesimal disk. 

3.1. Initial Conditions and Integration Times 

Each system begins with either iV iig = 5, 4, or 3 oli- 
garchs, together with the gas giants, Jupiter and Sat- 
urn. Initial conditions are the same as those described 
in H2.2.11 except that initial eccentricities and inclina- 
tions of oligarchs are such that ej — sinij = 10~ 2 . 
Initial inclinations are therefore ten times larger than 
the assumed inclination dispersion of disk planetesimals 
((sin 2 idisk) 1 ^ 2 = 10 -3 ). Initial eccentricities are the 
same as those that characterize phase 2b of the viscous 
stirring curves (Fig.QJi). Planetary orbits initially do not 
cross. 

Dynamical friction is exerted by a disk of constant sur- 
face density a which extends from a barycentric radius 
of 12.5 AU to 45 AU. The outer boundary coincides with 
the location of the classical Kuiper belt (C06). The in- 
ner boundary is less well motivated. It is chosen so that 
the oligarchs reside initially inside the disk while the gas 
giants do not. We explore values for a ranging from 0.4 
to 0.001 g/cm 2 . For reference, the initial surface density 
in oligarchs is S w 1.5 g/cm 2 . 
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In a subset of runs, we include 400 test particles repre- 
senting large KBOs. These feel the gravity of the plan- 
ets but do not feel dynamical friction from the disk (see 
£|2.ip . Initial semi-major axes of test particles range from 
a = 40 to 45 AU, and initial eccentricities and inclina- 
tions are such that e = sini = 10~ 2 . For all planets and 
KBOs, initial longitudes of ascending node, longitudes 
of pericenter, and mean anomalies are randomly chosen 
from uniform distributions between and 2ir. 

Settings for the MERCURY6 code are the same as 
those given in ^2.2. 1[ except for the duration of inte- 
gration. The integration automatically halts when there 
are a catastrophic number of ejections, i.e., when the 
only massive bodies remaining include the Sun, Jupiter, 
Saturn, and one oligarch (which may have collided with 
other oligarchs). In the absence of such an event, each 
system is integrated first to t = 2 x 10 7 yr. If the planets 
seem to have stabilized at that time — i.e., their eccentric- 
ities no longer grow — then we stop the integration and 
record the outcome as final. Otherwise, we repeat this 
test as necessary at 5 x 10 7 yr and 1 x 10 s yr. 

By t = 1 X 10 8 yr, most but not all realizations stabilize. 
Excluding systems that are stopped abruptly once only 
three planets remain, we find that < 10% of systems 
have undergone a close encounter (here defined to occur 
when the distance between any two planets is less than 
1 Hill radius) within the last 10 7 yr of the integration, 
for iVoiig = 5 and all values of a tested. For iV ug = 4 
and 3, the corresponding fractions are < 10% and < 
1%. Those realizations that do not stabilize by t = 1 x 
10 s yr are typically characterized by small values of a < 
0.006 g/cm 2 (i.e., relatively weak dynamical friction) and 
one oligarch remaining on an eccentric orbit that extends 
well past the outer edge of the disk. For these low values 
of a, the number of ejected oligarchs that we report will 
be underestimated, but not in a way that changes our 
qualitative conclusions. 

The initial conditions just summarized apply to all re- 
sults in the following two sections, § ^3.2rl3~3l An alter- 
native set of initial conditions, motivated by the findings 
in those sections, and results pertaining thereto are pre- 
sented in ^3.41 

3.2. Threshold Disk Surface Densities for Instability 

Oligarchs cross orbits when the disk surface density 
a is so low that dynamical friction cooling cannot com- 
pete with viscous stirring. GLS04 estimate the critical 
surface density for instability to be on the order of the 
surface density of oligarchs: <7 C rit ~ S (see also Chiang 
& Lithwick 2005 for a correction in the derivation of this 
result). How well does this criterion predict the onset of 
instability for our system of A i; g = 5 oligarchs? 

As noted in £ 12.21 the rates of eccentricity growth (vis- 
cous stirring) exhibited in our N-body integrations dif- 
fer from those estimated by GLS04. Specifically, our 
rates are slower, as evidenced by the period of slow-to- 
no growth of eccentricity (phase 2b) in Fig. QJi. Over- 
estimating the vigor of viscous stirring leads to overes- 
timates for tTcrit- We try to predict cr cr it ourselves by 
drawing from the numerical results of W2.2\ We observe 
in Fig. [IJi that after a time instable ~ 5 x 10 3 yr, ec- 
centricities surge rapidly to crossing values. 6 Therefore 

6 Chambers et al. (1996) provide fitting formulae for t uns table 
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Fig. 4. — Branching ratios for the outcome of the velocity in- 
stability, for systems that include dynamical friction with a disk 
of surface density a. For comparison, the initial surface density in 
oligarchs is E ss 1.5 g/cm 2 . Systems begin with Jupiter, Saturn, 
and either iV ii g = 5 (top panel), 4 (middle panel), or 3 (bottom 
panel) oligarchs. Different symbols denote the fraction of realiza- 
tions that result in 0, 1, 2, 3, or 4 ejections of oligarchs, whether 
collisionally merged or not. Few oligarchs collide in these simula- 
tions; the fraction of realizations that result in 1 collision is less 
than 20% for all values of a tested, and the fraction of realizations 
that result in > 1 collision is < 1%. To produce > 1 ejection with 
> 20% probability requires a < S/10. Standard error bars due to 
counting statistics are given for curves corresponding to iV i; g — 2 
ejections. 

for oligarchs to cross orbits, eccentricities must not be 
allowed to vanish by dynamical friction before £ uns tabie : 

^■vanish ^ ^unstable • 

(17) 

Using (|12| for i V anish, we find that (fl7|) translates into 
o~ < <7 C rit) where 

V2 m Q naef nit 
crit ~41n(l+A 2 ) m Gt unstahle [ ° J 

- 0.2 g/cm 2 - 0. IS. (19) 

The numerical evaluation takes a = 19.4 AU, emit = 0.01 
(the value appropriate to phase 2b of the viscous stirring 
curves in Fig. QJi; the median eccentricity does not rise 
above 0.01 for t < instable), and ln(l + A 2 ) = ln(l + 
A 2 nit ) = 0.02. 

Our semi-empirical estimate for cr cr i t finds support in 
Fig. UK, which documents, for all runs starting with 
A ng = 5 oligarchs, the frequency of incidence of out- 
comes ("branching ratios") as a function of o. Instability 
and the subsequent ejection of 1 and only 1 oligarch is 
the dominant outcome for a = 0.1 g/cm 2 « 0.072. For 
a > 0.3S, more than 90% of realizations produce no ejec- 
tion. Fig. [5] displays a sample simulation for a w 0.07E 
in which 1 oligarch escapes before the system stabilizes. 

Figs. HJa and 0fc supply branching ratios for A i; g = 4 
and 3. To produce the same number of ejections with 
smaller AT ij g (less viscous stirring) requires smaller a 
(less dynamical friction). For example, for N n g = 3, 
the ejection of 1 and only 1 oligarch is the dominant 

as functions of oligarch mass and orbital spacing, but only for the 
case without the gas giants Jupiter and Saturn. 
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Fig. 5. — Typical realization for a = 0.1 g/cm 2 ks 0.07S in which 
1 oligarch escapes before the system stabilizes. The top panel plots 
semi-major axis (solid curves), and periastron and apastron dis- 
tances (dotted curves) for each of the seven planets. Eccentricities 
(middle panel) and inclinations (bottom panel) are plotted using 
a mixed log-linear ordinate. Black curves refer to Jupiter, gray 
curves refer to Saturn, and remaining colored curves refer to the 
five oligarchs. 

outcome for a w 0.031], occurring in about 60% of real- 
izations. The corresponding a for N n g = 4 is 0.05S. 

3.3. Runs with Two Surviving Oligarchs 
(Solar-System-Like Outcomes) 

While (Tcrit ~ 0.1E roughly characterizes the onset of 
instability and the subsequent ejection of a single oli- 
garch, the disk surface density must be reduced below 
Ccrit to produce more than 1 ejection in a large fraction 
of runs. To generate an outcome reminiscent of our solar 
system starting with -/V ii g = 5 requires 3 ejections and 
the survival of 2 oligarchs. According to Fig. [4Jl, such an 
outcome occurs with a maximum probability of ^50% 
for a ~ 0.01 g/cm 2 w 0.007S. The probability exceeds 
20% for all values of a < 0.03S that we tested. 

Figs. and 0fc indicate that for N u g = 4 and 3, 
the values of a most likely to produce 4-planet sys- 
tems (Jupiter, Saturn, plus 2 oligarchs) are ~0.01S and 
~0.03E, respectively. The probabilities for generating 4- 
planet systems starting with N n g = 4 or 3 reach large 
maximum values of about 50%, and remain above 20% 
over a considerable range in a, up to ^0.07S in the case 
A^oiig = 3. 

The vast majority of the resultant 4-planet systems 
are correctly ordered; they contain, in order of increas- 
ing semi-major axis, Jupiter, Saturn, and 2 oligarchs. 
Moreover, in most of these systems, the surviving plan- 
ets have not experienced a collision. In the following 
sub-sections we further quantify the properties of these 
correctly ordered, collisionally unmodified, 4-planet sys- 
tems, comparing them to those of the solar system. We 
refer to the 2 surviving oligarchs in each of these systems 
as Uranus and Neptune. 

3.3.1. Final Semi-Major Axis Distributions 

Because packed oligarchies evolve chaotically, we can 
only meaningfully compute probability distributions for 
their final semi-major axes. Fig. [5] illustrates how these 
distributions depend on a, for realizations starting with 
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Fig. 6. — Cumulative distributions of final semi-major axes of sys- 
tems that contain, in order of increasing semi-major axis, Jupiter, 
Saturn, and two oligarchs (Uranus and Neptune), each with their 
original mass. These 4-planet systems result from integrations orig- 
inally containing 7 planets (N u g = 5). Each panel corresponds to 
simulations performed with the disk surface density a indicated. 
Vertical lines mark the current semi-major axes of the giant plan- 
ets in our solar system. Apart from a tendency for Saturn and 
Neptune to be pulled toward the inner (a = 12.5 AU) and outer 
(a = 45 AU) edges of the disk with increasing a (increasing dy- 
namical friction), the characteristic final semi-major axes do not 
depend on a. Saturn, Uranus, and Neptune have final orbits too 
large compared to their actual counterparts in the solar system. 

-Woiig — 5. Increasing a increases dynamical friction and 
therefore tends to pull Saturn and Neptune, whose or- 
bits lie near disk boundaries, into the disk. For example, 
if Saturn's orbital apocenter intersects the disk while its 
pericenter remains outside the disk, then dynamical fric- 
tion will circularize the orbit by raising the pericenter 
closer to apocenter. The kinks in the distribution func- 
tions for Saturn in Figs.[6K and[6b are located at a = 12.5 
AU, exactly at the inner disk edge. The kink vanishes in 
Fig. [Hh- For the simulations in Fig[5h, dynamical friction 
is strong enough to pull Saturn's orbit wholly into the 
disk at a > 12.5 AU. 

To improve statistics, we ignore these small artifacts of 
our disk boundary conditions and pool together N rea x = 
438 realizations, all of which start with 7V ii g = 5 and 
produce 4-planet systems, but have a variety of er's be- 
tween 0.001 and 0.1 g/cm 2 . From this pool we construct 
the distribution of final semi-major axes shown in Fig. [7^. 
Clearly, most realizations end with Saturn, Uranus, and 
Neptune on orbits too large compared to their actual 
counterparts in the solar system. Furthermore, Jupiter 
typically migrates inward as a consequence of ejecting 
several oligarchs outward. The excessively large orbits of 
Saturn, Uranus, and Neptune are a consequence of those 
planets having scattered oligarchs inward to Jupiter. 

Though they only comprise (given our assumed ini- 
tial conditions) a few percent of outcomes for a — 
0.02 g/cm 2 , some realizations better resemble the solar 
system insofar as Uranus and Neptune have final semi- 
major axes less than 30 AU. Fig.[8]showcases an example. 
Even for this simulation, however, Saturn's final orbit is 
3 AU larger than its actual one. 

Reducing the number of starting oligarchs significantly 
lessens the problem of excessive migration. Figs. [TJd and 
[7J; show final semi-major axis distributions corresponding 
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Fig. 7.— 



Cumulative distributions of final semi-major axes of 
correctly ordered, collision-free, 4-planet systems, shown for dif- 
ferent starting values of iV ij g . Each panel combines integrations 
performed with a variety of cr's between 0.001 and 0.1 g/cm 2 . Kinks 
in the distributions for Saturn and Neptune are artifacts of our as- 
sumed disk boundaries at a = 12.5 AU and 45 AU. Vertical lines 
mark the current semi-major axes of solar system giants. Decreas- 
ing 2V D ijg shrinks the final orbits for Saturn, Uranus, and Neptune. 
Final orbits resembling those of the current solar system are most 
easily produced for iV iig = 3. 
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Fig. 8. — Example realization for a = 0.02 g/cm 2 and N ^ s = 5 
in which 3 oligarchs (denoted by magenta, green, and cyan curves) 
escape before the system stabilizes. The remaining 2 oligarchs 
(blue and red) have final, nearly circular and co-planar orbits inside 
30 AU, similar to those of Uranus and Neptune. Such an outcome 
occurs with a frequency of a few percent for a = 0.02 g/cm 2 . 

to N \ ig = 4 and 3, respectively. Outcomes for -/V i ig = 3 
are most solar-system-like. In £13.41 we will experiment 
with initial conditions to demonstrate that a level of 
agreement comparable to that displayed in Fig. [7J; for 
N \i E = 3 can also be obtained for iVoiig = 4. 

3.3.2. Stirring of KBOs 

Fig. [9] describes how the test particles (KBOs), ini- 
tially distributed in a dynamically cold ring at a = 40- 
45 AU, are stirred by oligarchs, for the same iVoiig = 5 
simulation (Fig. [8|) which places Uranus and Neptune on 
final orbits inside 30 AU. Simulation data are collected 
at t = 5 x 10 7 yr. For comparison, Fig. [9] also plots data 
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Fig. 9. — Orbital elements of simulated KBOs (black crosses) at 
t = 5 X 10 7 yr, for the same simulation shown in Fig. \E\ Eccen- 
tricities and inclinations are displayed against either semi-major 
axis a or perihelion distance r p . For comparison we show orbital 
elements of known classical and scattered (non-resonant) KBOs 
(red inverted Y's) from C06. Orbital elements of surviving planets 
(Jupiter, Saturn, and the two oligarchs) are shown as blue circles. 
Simulated KBOs for this run do not exhibit the large eccen tric ities 
and inclinations characterizing actual KBOs, but see Figs. [TT1 and 
I14l for other examples. 

for actual KBOs that do not reside in any strong mean- 
motion resonance. These objects, taken from C06, com- 
prise both low-e classical and high-e scattered KBOs as 
classified by the Deep Ecliptic Survey (Elliot et al. 2005). 
The simulated test particles have their eccentricities and 
inclinations excited up to ~0.1, values that match actual 
classical KBOs. But the simulated particles fail to em- 
body the extreme degree of dynamical heating exhibited 
by scattered KBOs. Marauding oligarchs in this simula- 
tion stir planetesimals at 40-45 AU too briefly. 

Since simulations with iVoiig = 3 more efficiently gen- 
erate solar-system-like planetary spacings than do sim- 
ulations with iVoiig = 5, we can more thoroughly map 
out the possible extents to which KBOs are stirred for 
Aoiig = 3. Figs. [TU] and QT] document one simulation, 
representative of several percent of the solar-system-likc 
realizations generated using iVoiig = 3, in which KBOs 
are stirred considerably. Even here, however, the pro- 
portion of simulated KBOs that simultaneously attain 
inclinations i > 10° and perihelion distances r p > 35 
AU is less than observed. The proportion of simulated 
KBOs having eccentricities e > 0.3 also seems under- 
represented. 

3.4. More Compact Initial Conditions 

As described in ? ^3.2ff331 the simulations that begin 
with iVoiig = 4 or 5 oligarchs often do yield 4-planet sys- 
tems, but the orbital spacings of the resultant systems 
do not match those of the solar system. Multiple ejec- 
tions displace Jupiter too far inward and displace Sat- 
urn, Uranus, and Neptune too far outward. The prob- 
lem of excessive spreading is exacerbated by the need to 
have Neptune conclude its orbital evolution by migrating 
smoothly and slowly outward from a ~ 23 to 30 AU to 
produce the population of resonant KBOs (Murray-Clay 
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Fig. 10. — Example realization for N u B = 3 and a = 0.04 g/cm 2 
in which surviving planets have final orbits res emb ling those of the 
solar system (for an even closer match, see Fig. 1131 generated using 
a revised set of initial conditions). Moreover, the KBOs in this 
simulation, shown in Fig. 1111 are stirred considerably by oligarchs. 
Such an outcome represents several percent of runs generated using 
7V olig = 3. 



CT = 0.04g/cm a run=109 




10 100 
a (AU) 



5 10 50 
r„ (AU) 



Fig. 11. — Orbital elements of simulated KBOs at t = 5 X 10 7 yr, 
for the same simulation shown in Fig. 1101 In every panel except 
the one at the bottom left, simulated KBOs (black crosses) located 
to the left of observed KBOs (red inverted Y's) are unstable over 
the age of the solar system because of perturbations by the giant 
planets. Stable simulated KBOs (lying on top of and to the right 
of observed KBOs) exhibit large eccentricities and inclinations as a 
consequence of scattering by velocity-unstable oligarchs. Neverthe- 
less, the simulations fail to reproduce the most extreme of observed 
scattered KBOs having sin i > 0.2. 

& Chiang 2006, and references therein; but see Levison 
et al. 2006 for an alternative theory for the origin of reso- 
nant KBOs). This last constraint implies that our simu- 
lations should place Neptune on a final orbit near a ss 23 
AU. 

In this sub-section we attempt to remedy the problem 
of excessive spreading by adjusting our initial conditions. 
In anticipation of Jupiter's inward displacement, we lo- 
cate that planet initially at aj = 5.7 AU. In anticipation 



of Saturn's outward displacement, we set as = 8 AU ini- 
tially. The innermost oligarch is also shifted inwards, to 
ai = 12 AU. Initial semi-major axes for remaining oli- 
garchs are still given by eqn. (fTJ): 02 through 05 equal 
13.7, 15.5, 17.7, and 20.1 AU. Finally, so that all oli- 
garchs lie initially inside the disk, we extend the inner 
edge of the disk inward to 10 AU. All remaining param- 
eters remain unchanged from their values in Ej3.ll 

The distribution of final semi-major axes for resultant 
4-planet systems is given by Fig. 1121 constructed in simi- 
lar fashion to Fig. [7] The more compact initial configura- 
tion produces reasonably close matches to current orbital 
spacings in the solar system for N \ ig = 4. Results for 
A^oiig = 3 are also acceptable if we allow for the subse- 
quent outward migration of Neptune that is seemingly 
demanded by resonant KBOs (see first paragraph of this 
sub-section). The case iV ii g = 5 still suffers, however, 
from excessive spreading. 
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Fig. 12. — Same as Fig. [7] but gener ated using the more com- 
pact initial configuration described in £13.41 Runs starting with 
^olig = 5 oligarchs still produce final orbits for the 2 surviving 
oligarchs that are still too large compared to those of ice giants in 
our solar system. The case N a n g = 3 produces ice giant orbits that 
are too small. Nevertheless, Af i; g = 3 is acceptable if we allow 
for subsequent outward migration of the ice giants. Such outward 
migration seems necessary to explain the origin of resonant KBOs 
(C06; but see Levison et al. 2006 for an alternative theory). Re- 
sults for N n s = 4 are intermediate and produce orbital spacings 
most closely resembling those of the current solar system. Data for 
^olig = 5, 4, and 3 are generated with runs having a = 0.01, 0.04, 
and 0.1 g/cm 2 , respectively. 



Figs. [12] and Q3] sample one simulation using the re- 
vised compact configuration for N \i g = 4. We high- 
light this simulation because it reproduces solar system 
properties, insofar as (1) Uranus and Neptune have fi- 
nal semi-major axes less than 30 AU, and (2) the Kuiper 
belt at 40-45 AU is significantly stirred. Though out- 
come (1) is not infrequent — occurring in, e.g., 16 out of 
100 runs with a — 0.04 g/cm 2 and iVoiig = 4 — outcome 
(2) is less probable, characterizing only several percent of 
runs already culled to satisfy (1). Most runs that satisfy 
(1) stir KBOs to eccentricities and inclinations of just a 
few percent. By contrast, the simulation showcased in 
Fig. 1141 excites large eccentricities and inclinations simi- 
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lar to those sported by actual KBOs. Nevertheless, the 
most extreme of scattered KBOs, having perihelion dis- 
tances r p > 40 AU and inclinations i > 20°, are still 
under-represented. In short, our revised compact con- 
figuration stirs KBOs to about the same degree as our 
original configuration. 
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Fig. 13. — Similar to Fig. 1101 except generated for j V ij g = 4 
using the more compact initial configuration described in £|3.4I The 
simulation starts with a total of 6 planets (Jupiter, Saturn, and 4 
oligarchs) and ends with Jupiter, Saturn, and two oligarchs on 
orbits that closely match their actual counterparts in the solar 
system. 
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Fig. 14. — Orbital elements of simulated KBOs at t = 5 X 10 7 yr, 
for the same simulat ion shown in Fig. 1131 The same remarks given 
in the caption to Fig. 1 1 1 I concerning the stability of simulated KBOs 
apply here. Stable simulated KBOs, though stirred considerably 
by velocity-unstable oligarchs, fail to embody the large inclinations 
i > 20° exhibited by observed scattered KBOs. 



4. SUMMARY AND OUTLOOK 

In ij4.ll we answer the three questions posed in SJTJ In 
£ )4.21 we place our work in a broader context and mention 
some directions for future work. 



4.1. Answers to Questions Posed in J7] 

1. Of all our simulations that initially place iV iig = 5 
Neptune-mass planets between 15 and 25 AU, and 
that have disk surface densities a w 0.1 g/ cm 2 w 
0.07E (where E ~ 1.5g/cm 2 is the original sur- 
face density in oligarchs), 50% result in the ejection 
of a single oligarch (Fig. 0|. For runs that begin 
with Aoiig = 4 and 3 oligarchs, we achieve similar 
outcomes for cr/E ~ 0.05 and ~0.03, respectively. 
Jupiter is responsible for the vast majority of ejec- 
tions, which occur within ~10 7 yr. The likelihood 
of a single ejection remains as high as 20% if the 
above a's are increased by factors of 2-3. Roughly 
speaking then, we find that instability and ejection 
require er/E < 0.1. By comparison, GLS04 esti- 
mate that ct/E < 1 for instability. The difference 
arises partly because nearest-neighboring oligarchs 
in our simulations are separated by 5 Hill sphere 
radii, whereas their analysis of shear-dominated oli- 
garchy assumes the separation is closer to ~1 Hill 
sphere radius. Our choice of 5i?u is motivated by 
the half-width of an oligarch's annular feeding zone 
in a shear-dominated disk. This half-width spans 
2.5i?H (Greenberg et al. 1991). Because oligarchs 
separated by 5i?u viscously stir each other more 
slowly (Fig. [TJ phase 2) than do oligarchs separated 
by li?H, we find a threshold value for a lower than 
what GLS04 estimate. 

2. For certain choices of a and initial semi-major axes, 
systems starting with N u g = 3 or 4 oligarchs fre- 
quently end with 2 surviving oligarchs on nearly 
circular and co-planar orbits inside 30 AU (Figs. [7] 
and [T2"|) . For example, of all runs that (a) use 
our revised set of initial semi- major axes fi )3.4p . 
(b) begin with A i; g = 4 oligarchs, and (c) have 
a = 0.04 g/cm 2 w 0.02E, 44% end with solar- 
system-like configurations in which the outermost 
surviving oligarch orbits inside 30 AU. This per- 
centage decreases with increasing A i; g . This is 
because surviving oligarchs spread outward, well 
beyond the current orbit of Neptune, as they scat- 
ter more oligarchs inward for eventual ejection by 
Jupiter. To eject efficiently more than one oligarch 
requires that a be reduced considerably below the 
previously mentioned threshold of ^O.IE. For ex- 
ample, we find that for A ii g = 4 and our origi- 
nal set of initial semi- major axes ( ^2.2. ip . setting 
a ~ 0.02 g/cm 2 ~ 0.01E maximizes the likelihood 
of 2 ejections at ~50%. 

3. In a small fraction of runs that successfully place 
Jupiter, Saturn, and 2 oligarchs on solar-system- 
like orbits inside 30 AU, test particles (KBOs) lo- 
cated initially in a dynamically cold ring at 40- 
45 AU have their eccentricities and inclinations 
considerably excited by velocity- unstable oligarchs. 
We observe maximum eccentricities of ~0.8 and 
maximum inclinations of ~20° (Figs. ITTI and fT4|). 
In runs characterized by the greatest degrees of ex- 
citation, orbits of simulated KBOs resemble those 
of observed classical KBOs and some observed scat- 
tered KBOs. However, no run reproduces the large 
proportion of observed scattered KBOs having in- 
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clinations > 20°. There may also be a problem 
in generating enough KBOs having eccentricities 
> 0.3, given the observational selection bias against 
finding such objects. 

These results are complementary to those of Levison 
& Morbidelli (2007), who concentrate on the limit a ^> 
E and find that they cannot produce solar-system-like 
outcomes. In comparison, we study the case a < £ and 
find positive results. 

4.2. Commentary 

4.2.1. The Compactness of Our Preferred Initial Conditions 

We find that a shear-dominated oligarchy can readily 
produce a solar-system-like outcome if it contains just 
a few excess oligarchs — about 1 or 2 extra — and if the 
oligarchs initially reside inside 20 AU. We are driven to 
these parameters because to scatter excess oligarchs in- 
ward (toward Jupiter for eventual ejection), surviving 
oligarchs must be scattered outward. The right amount 
of outward spreading is achieved for suitably compact 
initial configurations and not too many ejections. 

Our favored initial conditions are about as compact 
as those of Tsiganis et al. (2005), who place their out- 
ermost ice giant initially at 17 AU. By comparison, in 
our revised set of initial conditions for A ii g = 4, the 
outermost oligarch is located at 17.7 AU. But we stress 
that our study differs from theirs in that we base our 
initial conditions on considerations of shear-dominated 
oligarchic accretion. An ice giant cannot form at 17 
AU within the gas photoevaporation time of a few x 
10 7 yr without strong gravitational focussing (Sjl] Levi- 
son & Stewart 2001; GLS04). This need for gravitational 
focussing can be met by a highly dissipative disk of plan- 
etesimals (GLS04; Rafikov 2004). It is this disk, and the 
multiple (> 2) ice giants that it spawns, that we have 
modeled. 

Why such a disk would not form Neptune-mass oli- 
garchs outside 20 AU is an open question. For some 
ideas on what limits the sizes of planetary systems, see 
Youdin & Shu (2002) and Youdin & Chiang (2004). 

4.2.2. Reducing the Disk Surface Density ("Clean-Up" and 

Migration) 

Velocity instability and the ejection of a single oligarch 
require er/E < 0.1. Ejecting more than 1 oligarch re- 
quires still lower values of c/E ~ 0.01. How can a reach 
such low values? Accretion and/or ejection of planetes- 
imals by oligarchs are natural possibilities. The rate at 
which a decreases, compared to the rate at which oli- 
garchs stir each other, determines whether more than 1 
oligarch escapes. If the rate of depletion of a is suffi- 
ciently slow, then after the first oligarch escapes, sur- 
viving planets may occupy orbits so spread apart that 
they remain stable even as a decreases further. On the 
other hand, if the rate of depletion is fast, then condi- 
tions required to eject more than 1 oligarch can be met. 
We leave investigation of the time dependence of a for 
future work. 

Given an initial surface density in oligarchs of E ~ 
1 g/cm 2 , the disk surface densities relevant for instability 
and ejection range from a « 0.1 to 0.01 g/cm 2 . These cr's 
are still higher than surface densities characterizing the 



Kuiper belt today, which are on the order of 0.001 g/cm 2 
(integrated over all KBO sizes). "Cleaning up" the disk 
mass remains an unsolved problem (GLS04). Again, the 
mass could either be accreted or ejected by surviving 
planets. Planetesimal ejection drives planetary migra- 
tion. For Neptune to expand its orbit from ^23 to 30 
AU, as seemingly demanded by the large observed popu- 
lation of resonant KBOs (C06), the planet must scatter 
at least ~7/30 ~ 25% of its own mass in planetesimals, 
or about AMq. Spreading such a mass over a disk of 
radius 30 AU yields a surface density on the order of 
0.03 g/cm 2 , which falls within the range of cr's that we 
find characterize instability. Clean-up and migration go 
hand-in-hand. 

The disk masses that we found to be relevant for insta- 
bility are smaller than the masses in the planets. This 
raises concern about the validity of our approximation 
that disk properties remain fixed throughout the simula- 
tion. One way of testing this approximation is to com- 
pare the change in the total angular momentum of all 
planets (brought about by dynamical friction) to the an- 
gular momentum available in the disk. If the former 
were much larger than the latter, then our neglect of 
back-reaction upon the disk would be a poor assump- 
tion. For the simulations displayed in Figs. [10] and [T3l 
we find that the change in the z-component of angu- 
lar momentum of all planets (including ejected ones) is 
nearly identical to the angular momentum available in 
the disk, indicating that our assumption of a fixed disk 
may be only marginally valid. (An analogous test for the 
energy would be inconclusive since the disk is supposed 
to be a sink of energy by virtue of dissipative collisions). 
For ideas on how to treat planet-disk interactions self- 
consistently, see Lithwick & Chiang (2007) and Levison 
& Morbidelli (2007). 

4.2.3. Disk Optical Depths and Comparison to Debris Disks 

A planetesimal disk of surface density a has 
a geometric vertical optical depth r p ~ 4 x 
10~ 6 [cr/(0.01g/cm 2 )](10m/p), where p is the assumed 
planetesimal radius. Collisions between planetesimals, 
which occur over a timescale ~l/(f2r p ) ~ 10 7 yr at 
30 AU, generate smaller dust particles whose opti- 
cal depth is orders of magnitude higher. For exam- 
ple, if the dust size distribution obeys a Dohnanyi 
(1969) spectrum, then the geometric, vertical optical 
depth in s-sized grains would be ~ T p (p/s) 1 ^ 2 ~ 
0.01[o-/(0.01g/cm 2 )](10m/p) 1 / 2 (/im/s) 1 / 2 . This is com- 
parable to the vertical optical depths of some of the 
brightest extra-solar debris disks observed, e.g., j3 Pic 
(Artymowicz 1997), HR 4796A (Li & Lunine 2003), and 
AU Mic (Strubbe & Chiang 2006), systems that are all 
~10 7 yr old. The observed paucity of stars with optically 
thicker disks implies that large populations of planetes- 
imals having sizes p < 10 m at stcllocentric distances of 
~30 AU cannot be maintained for longer than ^10 7 yr. 
The surface density in such collisional objects must be 
reduced by at least 2 orders of magnitude below planet- 
forming values of ~lg/cm 2 within this timescale. In 
other words, planetesimals that are both collisional and 
planet-forming, like the kind espoused by GLS04, must 
be cleaned up fairly quickly. 



14 



Ford & Chiang 



4.2.4. Instability in Our Solar System and Others 

We have demonstrated that more than 2 ice giants may 
once have orbited the Sun. The current architecture of 
the outer solar system may well have resulted from a 
prior era of dynamical instability during which Uranus, 
Neptune, and 1 or 2 other ice giants crossed paths. 

We expect similar instabilities to afflict all nascent 
planetary systems. Perhaps planet-planet instabilities 
are reflected in the large orbital eccentricities exhib- 
ited by most extra-solar gas giants (Marzari & Weiden- 
schilling 2002; Ford et al. 2003; but for an alternative 
view see Goldreich & Sari 2003). The case of Upsilon 
Andromedae fits this picture (Ford et al. 2005). The dif- 
ference between our solar system and systems like Ups 
And might be the surface density of the parent disk at 
the time of the last planet-planet scattering (Ford 2006). 
The time of last scattering will vary widely because of 
the chaotic nature of multi-planet systems. In the case of 
the solar system, the disk surface density must have been 
large enough at the time of last scattering for dynami- 
cal friction to damp the eccentricities and inclinations of 
surviving planets back down. 

Viscous stirring rates vary with the semi-major axis 
separation between oligarchs. At least in the case with- 
out gas giants, the time for oligarchs to undergo close 
encounters increases by several orders of magnitude as 
their semi-major axis spacing is increased from 3 to 7 mu- 
tual Hill radii (Chambers et al. 1996). The disk surface 
density required for instability depends directly on this 
time, i.e., cr cr it <x l/£ U nstabie in our eqn. (fTg]) . Whether 
^unstable varies as strongly with oligarch spacing when 
perturbations by Jupiter and Saturn are included is not 
known, but preliminary experiments by us suggest that 
it does not. When we change the oligarch spacing from 
our standard 5 Hill radii to 3 Hill radii in runs that in- 
clude Jupiter, Saturn, and N \ ig = 5 oligarchs, we find 
that the probability of 1 ejection still peaks at ^50% for 
cr/E w 0.1 (accounting for the factor of 2 increase in S 
due to the shorter spacing). 

4.2.5. Evidence for the Velocity Instability in the Kuiper 

Belt 

Did velocity-unstable ice giants excite the large eccen- 
tricities and inclinations of the scattered Kuiper belt, as 
proposed by C06? Our provisional answer is no, as we 



were unable to reproduce the large inclinations of scat- 
tered KBOs. For runs with disk surface densities down 
to ^0.01 g/cm 2 , oligarchs spend too little time, less than 
~10 7 yr, passing through the Kuiper belt. Moreover, oli- 
garchs in our simulations have orbital inclinations that 
rarely exceed 10°. 

To remedy the situation, we might appeal to still lower 
disk surface densities, on the order of ^0.001 g/cm 2 , 
for which dynamical friction cooling times for embed- 
ded planets would be as long as ~10 8 yr. We found this 
region of parameter space difficult to explore. Out of 
100 simulations starting with (a) iV ug = 4, (b) our com- 
pact set of initial conditions, and (c) a — 0.002 g/cm 2 , 
only 6 yielded systems each with 2 surviving oligarchs at 
the end of the integrations at t = 10 8 yr. Unfortunately, 
5 of these 6 systems had not stabilized, and it was un- 
clear whether more oligarchs would be ejected were the 
integrations to continue. Furthermore, over these long 
timescales, effects resulting from time variations in a (see 
§4.2 . 2() might be expected to be important, and we have 
not modeled these. 

Cooling times for planets and, by extension, KBO 
heating times might also be prolonged in more realistic 
treatments of dynamical friction that incorporate non- 
axisymmetries and the clearing of gaps in the disk. Ways 
of numerically simulating the response of planetesimal 
disks to planets can be found in Levison & Morbidelli 
(2007) and Lithwick & Chiang (2007). 



We thank Ben Collins, Yoram Lithwick, Ruth Murray- 
Clay, Re'em Sari, and an anonymous referee for helpful 
exchanges. We also thank Harold Levison and Alessan- 
dro Morbidelli for generously sharing some of their ideas 
about planet-disk interactions and suggesting that we 
examine the angular momentum and energy budgets of 
our simulations. Support for E.B.F. was provided by a 
Miller Research Fellowship and by NASA through Hub- 
ble Fellowship grant HST-HF-01195.01A awarded by the 
Space Telescope Science Institute, which is operated by 
the Association of Universities for Research in Astron- 
omy, Inc., for NASA, under contract NAS 5-26555. E.C. 
acknowledges grants from the National Science Founda- 
tion (AST-0507805), NASA (JPL-1264475), and the Al- 
fred P. Sloan Foundation. 



REFERENCES 



Artymowicz, P. 1997, Annual Review of Earth and Planetary 
Sciences, 25, 175 

Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton: 
Princeton University Press), 424 
Butler, R.P., et al. 2006, ApJ, 646, 505 
Chambers, J.E. 2006, Icarus, 180, 496 
Chambers, J.E. 1999, MNRAS, 304, 793 

Chambers, J.E., Wetherill, G.W., & Boss, A. P. 1996, Icarus, 119, 
261 

Chiang, E.I., & Lithwick, Y. 2005, ApJ, 628, 520 

Chiang, E.I., et al. 2006, in Protostars and Planets V, eds. B. 

Rcipurt h, D. Jew itt, & K. Keil (Tucson: Univ. Arizona Press), in 

press l |astro-ph/0601654} (C06) 

Collins, B.F., & Sari, R. 2006, A J, 132, 1316 

Collins, B.F., Schlichting, H.E., & Sari, R. 2007, AJ, in press 

Dohnanyi, J.W. 1969, JGR, 74, 2531 



Elliot, J.L., et al. 2005, AJ, 129, 1117 

Ford, E.B. 2006, in New Horizons in Astronomy: Frank N. Bash 
Symposium, ASP Conference Series Vol. 352 (San Francisco: 
Astronomical Society of the Pacific), 15 

Ford, E.B., Lystad, V., & Rasio, F.A. 2005, Nature, 434, 873 
Ford, E.B., Rasio, F.A., & Yu, K. 2003, in Scientific Frontiers in 
Research on Extrasolar Planets, ASP Conference Series Vol. 294, 
eds. D. Deming & S. Seager (San Francisco: Astronomical Society 
of the Pacific), 181 

Goldreich, P., Lithwick, Y., & Sari, R. 2004, ARA&A, 42, 549 
(GLS04) 

Goldreich, P., & Sari, R. 2003, ApJ, 585, 1024 
Goldreich, P., & Tremaine, S. 1982, ARA&A, 20, 249 
Greenberg, R., et al. 1991, Icarus, 94, 98 
Ida, S., & Makino, J. 1993, Icarus, 106, 210 
Kenyon, S.J., & Luu, J.X. 1999, 118, 1101 



Ice Giants 



Levison, H.F., & Stewart, G.R. 2001, Icarus, 153, 224 

Levison, H.F., et al. 2006, in Protostars and Planets V, eds. B. 

Reipurth, D. Jewitt, & K. Keil (Tucson: Univ. Arizona Press), in 

press 

Levison, H.F., & Morbidelli, A. 2007, Icarus, in press 
Li, A., & Lunine, J. I. 2003, ApJ, 590, 368 

Lissauer, J.J., Pollack, J.B., Wetherill, G.W., & Stevenson, D.J. 
1995, in Neptune and Triton, ed. D. Cruikshank (Tucson: Univ. 
Arizona Press), 37 

Lithwick, Y., & Chiang, E. 2007, ApJ, in press (astro-ph:0607241) 
Marzari, F., & Weidenschilling, S.J. 2002, Icarus, 156, 570 
Matsuyama, I., Johnstone, D., & Hartmann, L. 2003, ApJ, 582, 
893 

Murray-Clay, R.A., & Chiang, E.I. 2006, ApJ, 651, 1194 
Rafikov, R.R. 2004, AJ, 128, 1348 

Shu, F.H., Johnstone, D., Hollenbach, D. 1993, Icarus, 106, 92 

Stewart, G.R., & Ida, S. 2000, Icarus, 143, 28 

Strubbe, L.E., & Chiang, E.I. 2006, ApJ, 648, 652 

Thommcs, E.W., Duncan, M.J., & Levison, H.F. 1999, Nature, 

402, 635 

Thommes, E.W., Duncan, M.J., & Levison, H.F. 2002, AJ, 123, 
2862 

Tsiganis, K., et al. 2005, Nature, 435, 459 
Wisdom, J., & Holman, M. 1991, AJ, 102, 1528 
Youdin, A.N., & Chiang, E.I. 2004, ApJ, 601, 1109 
Youdin, A.N., & Shu, F.H. 2002, ApJ, 580, 494 



